__global__ void Interval3DDDA2(ABSSURF_SHARE_INFO asi, float* RayDir, uint* Depth3DDDA, uint* Interval, int maxIntel, int* PatchIndex, int nPatch, uint pixelBits, uint patchBits, unsigned char* surface, size_t pitch, int onlyClimb, float4 gPlane, int Shadow, unsigned char* surfNorm)
{
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=cRenderParam.renderWidth || y>=cRenderParam.renderHeight)			return;
	uint pixelI=x+y*cRenderParam.renderWidth;

	float* pixel = (float*)(surface + y*pitch) + 4*x;
	float dir[3];
	float eye[3];
	if(Shadow==0)
	{
		pixel[0]=pixel[1]=pixel[2]=1.0f;
		eye[0] = cRenderParam.eye[0];
		eye[1] = cRenderParam.eye[1];
		eye[2] = cRenderParam.eye[2];   //useful for tranparent illustration

		float screen[3] ={2*(x+0.5f)/cRenderParam.renderWidth-1.0f, 1.0f-2*(y+0.5f)/cRenderParam.renderHeight, 0.0f};
		float trans[4];
		BackProjGPU(screen, trans);
		dir[0] = trans[0]/trans[3]-eye[0];
		dir[1] = trans[1]/trans[3]-eye[1];
		dir[2] = trans[2]/trans[3]-eye[2];


		float t = -(gPlane.x*eye[0]+gPlane.y*eye[1]+gPlane.z*eye[2]+gPlane.w)/(gPlane.x*dir[0]+gPlane.y*dir[1]+gPlane.z*dir[2]);
		float pt[3]={eye[0]+t*dir[0],eye[1]+t*dir[1],eye[2]+t*dir[2]};

		if(fabs(pt[0])<GROUND_CORNER&&fabs(pt[2])<GROUND_CORNER)
		{
			//float xxf=(pt[0]+GROUND_CORNER)*DELTA_DIS,yyf=(pt[2]+GROUND_CORNER)*DELTA_DIS;
			AAChessBoard(screen, gPlane, eye, pt[0], pt[2], pixel);

			eye[0] = cRenderParam.light[0].pos[0];
			eye[1] = cRenderParam.light[0].pos[1];
			eye[2] = cRenderParam.light[0].pos[2];   //useful for tranparent illustration
			dir[0] = pt[0]-eye[0];
			dir[1] = pt[1]-eye[1];
			dir[2] = pt[2]-eye[2];
		}
		else
		{
			Depth3DDDA[pixelI] = 0;
			return;
		}
	}
	if(Shadow==1)
	{
		eye[0] = cRenderParam.eye[0];
		eye[1] = cRenderParam.eye[1];
		eye[2] = cRenderParam.eye[2];   //useful for tranparent illustration
		float screen[3] ={2*(x+0.5f)/cRenderParam.renderWidth-1.0f, 1.0f-2*(y+0.5f)/cRenderParam.renderHeight, 0.0f};
		float trans[4];
		BackProjGPU(screen, trans);
		dir[0] = trans[0]/trans[3]-eye[0];
		dir[1] = trans[1]/trans[3]-eye[1];
		dir[2] = trans[2]/trans[3]-eye[2];
	}
#if SELF_SHADOW
	if(Shadow==2)
	{
		if(fabs(pixel[3]-1.0)<1.0e-6 || fabs(pixel[3]-RENDER_AA)<1.0e-6)
		{
			float* pixelN = (float*)(surfNorm + y*pitch) + 4*x;
			eye[0] = cRenderParam.light[0].pos[0];
			eye[1] = cRenderParam.light[0].pos[1];
			eye[2] = cRenderParam.light[0].pos[2];
			dir[0] = pixelN[0]-eye[0];
			dir[1] = pixelN[1]-eye[1];
			dir[2] = pixelN[2]-eye[2];
		}
		else
		{
			Depth3DDDA[pixelI] = 0;
			return;
		}
	}
#endif

	//if(onlyClimb==0)
	uint pixelOffset = RAY_DIR_SIZE*pixelI;
	RayDir[pixelOffset  ] = dir[0];
	RayDir[pixelOffset+1] = dir[1];
	RayDir[pixelOffset+2] = dir[2];




	float tmin, tmax;
	int firstSeg[3] = {0, 0, 0};
	if(RayBoxIntersect(&tmin, &tmax, eye, dir, asi.box0, asi.boxSize, firstSeg, asi.seg)==0 
		|| ((Shadow==1)&&tmax<1) 
		|| (Shadow==0 &&(tmin<=0||tmax>=1.0)) 
		|| (Shadow==2 &&(tmin<=0))  )				//in selfshadow tmax>1.0
	{
		//if(onlyClimb==0)
		{
			//pixel[3] = 0.0f; 
			Depth3DDDA[pixelI] = 0; 	
		}
		return;
	}

	
	int patInd[3];				//patch index in traversal
	float pmin[3];
	tmin = tmin+RAY_ADD_EPS;
	pmin[0] = eye[0]+tmin*dir[0];
	pmin[1] = eye[1]+tmin*dir[1];
	pmin[2] = eye[2]+tmin*dir[2];
	patInd[0] = floorf((pmin[0]-asi.box0[0])/asi.boxSize[0]);
	patInd[1] = floorf((pmin[1]-asi.box0[1])/asi.boxSize[1]);
	patInd[2] = floorf((pmin[2]-asi.box0[2])/asi.boxSize[2]);

	int Di=0, Dj=1, mainD=2;			//main traversal direction is stored in the 3rd of mainD.
	int incDi, incDj, incMainD;			//increment diretion of every one.
	float slopeDi, slopeDj;				//the other 2 slope relative to the main dirction.
	float diffDi, diffDj;				//the other 2 diff
	float minDiff;

	GetMainDir(dir, &Di, &Dj, &mainD);
	float hi,hj,hmainD;
	if(dir[Di]<0)
	{
		incDi=-1;
		hi = pmin[Di]-asi.box0[Di]-patInd[Di]*asi.boxSize[Di];
	}
	else
	{
		incDi=1;
		hi = asi.box0[Di]+(patInd[Di]+1)*asi.boxSize[Di]-pmin[Di];
	}
	if(dir[Dj]<0)
	{
		incDj=-1;
		hj = pmin[Dj]-asi.box0[Dj]-patInd[Dj]*asi.boxSize[Dj];
	}
	else
	{
		incDj=1;
		hj = asi.box0[Dj]+(patInd[Dj]+1)*asi.boxSize[Dj]-pmin[Dj];
	}
	if(dir[mainD]<0)
	{
		incMainD=-1;
		hmainD = (pmin[mainD]-asi.box0[mainD])/asi.boxSize[mainD]-patInd[mainD];
	}
	else
	{
		incMainD=1;
		hmainD = (asi.box0[mainD]-pmin[mainD])/asi.boxSize[mainD]+(patInd[mainD]+1);
	}
	slopeDi = fabsf(dir[mainD]/dir[Di]);
	slopeDj = fabsf(dir[mainD]/dir[Dj]);
	diffDi = hi*slopeDi/asi.boxSize[mainD];
	diffDj = hj*slopeDj/asi.boxSize[mainD];
	minDiff = min(diffDi,diffDj);


	//trav3DDDA();
	int intelN=0, patchI;
	uint intervalOff=maxIntel*pixelI;
	while(intelN<maxIntel)
	{
		patchI = patInd[0]*asi.segShift[0]+patInd[1]*asi.segShift[1]+patInd[2]*asi.segShift[2];
		patchI = PatchIndex[patchI];
		if(patchI>=0&&patchI<nPatch)
		{
			//if(onlyClimb==0)
				Interval[intervalOff++] = (((x<<pixelBits)|y)<<patchBits)|(uint)(patchI);
			intelN++;
		}

		if(hmainD<minDiff)
		{
			patInd[mainD] += incMainD;
			hmainD += 1;
			if(patInd[mainD]>=asi.seg[mainD]||patInd[mainD]<0)	break;
		}
		else if(diffDi<diffDj)
		{
			patInd[Di] += incDi;
			diffDi += slopeDi;
			if(patInd[Di]>=asi.seg[Di]||patInd[Di]<0)			break;
		}
		else
		{
			patInd[Dj] += incDj;
			diffDj += slopeDj;
			if(patInd[Dj]>=asi.seg[Dj]||patInd[Dj]<0)			break;
		}
		minDiff = min(diffDi,diffDj);
	}

	
	Depth3DDDA[pixelI] = intelN;

}

__global__ void PeelSolver2(unsigned char* surfNorm, unsigned char* surface, size_t pitch, int degSum, uint* IntervalPad, float* RayDir, uint sortPadNum, ABSSURF_SHARE_INFO asi, uint pixelBits, uint patchBits, int* PatchIndexInverse, float* FuncData, uint* Unfinished, uint depthI, int Shadow)
{
	int thid = threadIdx.x;
	int thidA = thid+__mul24(blockIdx.x, blockDim.x);

	__shared__ uint patchInfo[PIXELTHREADS];
	patchInfo[thid] = IntervalPad[thidA];
	__syncthreads();

	int patchIndex = (int)( (patchInfo[0]<<(pixelBits+pixelBits))>>(pixelBits+pixelBits) );
	int segI = PatchIndexInverse[patchIndex];
	int segInd[3], segIndd[3];
	segInd[0] = segI/asi.segShift[0];
	segInd[1] = segI%asi.segShift[0]/asi.segShift[1];
	segInd[2] = segI%asi.segShift[1];
	segIndd[0] = segInd[0]+1;
	segIndd[1] = segInd[1]+1;
	segIndd[2] = segInd[2]+1;

	extern __shared__ float func[];
	int funcSize = (asi.deg[0]+1)*(asi.deg[1]+1)*(asi.deg[2]+1);

	if(thid<funcSize)
		func[thid] = FuncData[patchIndex*funcSize+thid];
	__syncthreads();

	if(patchInfo[thid]==0&&thid!=0)				{Unfinished[thidA]=0;return;}

	uint x = (patchInfo[thid]>>(patchBits+pixelBits));
	uint y = (patchInfo[thid]<<pixelBits)>>(patchBits+pixelBits);
	uint pixelI = x + y*cRenderParam.renderWidth;


	int indOff;
	indOff= RAY_DIR_SIZE*pixelI;
	float dir[3] = {RayDir[indOff], RayDir[indOff+1], RayDir[indOff+2]};
	float tmin, tmax;
	RayBoxIntersect(&tmin, &tmax, cRenderParam.light[0].pos, dir, asi.box0, asi.boxSize, segInd, segIndd);	
	//if((tmin<=0||tmax>=1))										{Unfinished[thidA]=0;return;}
#if SELF_SHADOW
	tmax = min(1.0,tmax);
	if(tmax<=tmin)		{Unfinished[thidA]=0;return;}
#endif


	float fValue[MAX_DEG_INTERPOLATION+1], coeff[MAX_DEG_INTERPOLATION+1];

	//********* code surround by /**/ is chebyshev interpolation ****************//
	//remember: change the interpolation here will also need change the InvMatrix;
	float In[3]={(cRenderParam.light[0].pos[0]+tmin*dir[0]-asi.box0[0])/asi.boxSize[0]-segInd[0],
		(cRenderParam.light[0].pos[1]+tmin*dir[1]-asi.box0[1])/asi.boxSize[1]-segInd[1], 
		(cRenderParam.light[0].pos[2]+tmin*dir[2]-asi.box0[2])/asi.boxSize[2]-segInd[2]};

	float Out[3];
	Out[0] = In[0]+(tmax-tmin)*dir[0]/asi.boxSize[0];
	Out[1] = In[1]+(tmax-tmin)*dir[1]/asi.boxSize[1];
	Out[2] = In[2]+(tmax-tmin)*dir[2]/asi.boxSize[2];

	float xyz[3];
	float cheby;
	for(int i=0; i<=degSum; i++)
	{
		cheby = (1-(__cosf(PI*(2*i+1)/(2*(degSum+1)))))/2;
		xyz[0] = (1-cheby)*In[0]+cheby*Out[0];
		xyz[1] = (1-cheby)*In[1]+cheby*Out[1];
		xyz[2] = (1-cheby)*In[2]+cheby*Out[2];
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, fValue+i);
	}
	//********* code surround by /**/ is chebyshev interpolation ****************//

	indOff = 0;
	for(int i=0; i<=degSum; i++)
	{
		coeff[i] = 0.0;
		for(int j=0; j<=degSum; j++)
		{
			coeff[i] += fValue[j]*cInvMatrix[indOff++];
		}
	}
	int indC;
	for(indC=0; indC<degSum; indC++)
	{
		if(coeff[indC]*coeff[indC+1]<0)	break;
	}
	if(indC==degSum)	{Unfinished[thidA]=1;return;}
	else
	{
#if SELF_SHADOW
		if(Shadow==2)
		{
#if _KI_BEZ_
		float firstRoot = KnotInsertBez(coeff, degSum);
#else
		float firstRoot = KnotInsert(coeff, degSum);
#endif
		if(firstRoot<0)		{Unfinished[thidA]=1;return;}

		firstRoot = tmin*(1-firstRoot) + tmax*firstRoot;
		float tPos = (*((float*)(surfNorm + y*pitch) + 4*x) - cRenderParam.light[0].pos[0])/dir[0];
		if(firstRoot<tPos*0.9)
		{
			float* pixel = (float*)(surface + y*pitch) + 4*x;
			pixel[0] = 0.8*pixel[0];
			pixel[1] = 0.8*pixel[1];
			pixel[2] = 0.8*pixel[2];	
			pixel[3] = SHADOW_RAY;						
		}Unfinished[thidA]=0;return;
		}
#endif
		float* pixel = (float*)(surface + y*pitch) + 4*x;
		pixel[0] = 0.4*pixel[0];
		pixel[1] = 0.4*pixel[1];
		pixel[2] = 0.4*pixel[2];	
		pixel[3] = SHADOW_RAY;
		
		Unfinished[thidA]=0;return;
	}	
}


extern "C" void KIABSPeelShadow(RENDER_PARAM* renderParam, unsigned char* surfNorm, unsigned char* surface, size_t pitch, bool* reNormal)
{
	cudaMemcpyToSymbol(cRenderParam, renderParam, sizeof(RENDER_PARAM));
	CUDA_SAFE_CALL( cudaMemset(surfNorm, 0, pitch*height) );
	CUDA_SAFE_CALL( cudaMemset(surface, 0, pitch*height) );


	//Computer shadow 3DDDA
	dim3 Db = dim3(16,16);
	dim3 Dg = dim3( width/Db.x+1, height/Db.y+1 );
	if(width %Db.x==0)	Dg.x-=1;
	if(height%Db.y==0)  Dg.y-=1;	
	Interval3DDDA2<<<Dg, Db>>>(asi, dRayDir, dDepth3DDDA, dInterval, maxInterval, dPatchIndex, nPatch, pixelBits, patchBits, (unsigned char*)surface, pitch, (int)onlyClimb, groundPlane, 0, surfNorm);
	
	uint scanDepth;
	uint depthI = 0;
	uint unfinishedN = pixelN;
	uint sortPadNum=0;
	dim3 grid = dim3(1,1,1);

	//printf("____________The shadow test_______________________________Begin____\n");
	CUDA_SAFE_CALL( cudaMemset(dIntervalCompact, 0, pixelN*sizeof(uint)) );
	while(true&&depthI<renderParam->nIter)
	{		
		grid.x = unfinishedN/SCAN_THREADN+1;
		if(unfinishedN%SCAN_THREADN==0)	grid.x-=1;
		CUDA_SAFE_CALL( cudaMemset(dScanDepth, 0, sizeof(uint)) );
		PeelDepth<<<grid, SCAN_THREADN, 1024>>>(unfinishedN, dDepth3DDDA, dInterval, maxInterval, dIntervalCompact, dScanDepth, depthI, pixelBits, patchBits);
		depthI++;

		CUDA_SAFE_CALL(cudaMemcpy(&scanDepth, dScanDepth, sizeof(uint), cudaMemcpyDeviceToHost) );
#if INFO_OUTPUT
		printf("%u step peel %u pixels.\n", depthI, scanDepth);
#endif	

		if(scanDepth==0)	break;	

		//Pad the interval to PIXELTHREADS-aligned
		dIntervalPad = dIntervalPeel;
		sortPadAll(scanDepth, &sortPadNum,  pixelN*3/2*sizeof(uint), dIntervalCompact, dIntervalPad, dIntervalPad, hIntervalSort, hIntervalPad, PIXELTHREADS);
#if INFO_OUTPUT
		printf("actual pad: %u, pad allocate: %u\n", sortPadNum, pixelN*3/2);
#endif
		//dumpBuffer(dIntervalPad, sortPadNum);

		if(sortPadNum==0)	break;
		if((int)(sortPadNum)>pixelN*maxAvgInterval)
		{
			printf("sortPad number > interval compact!\n");
			return;
		}

		//Hit test
		grid.x = sortPadNum / PIXELTHREADS + 1;
		if(sortPadNum%PIXELTHREADS == 0)		grid.x -= 1;
		CUDA_SAFE_CALL( cudaMemset(dUnfinished, 0, pixelN*3/2*sizeof(uint)) );
		PeelSolver2<<<grid, PIXELTHREADS, 1024>>>((unsigned char*)surfNorm, (unsigned char*)surface, pitch, degSum, dIntervalPad, dRayDir, sortPadNum, asi, pixelBits, patchBits, dPatchIndexInverse, dFuncData, dUnfinished, depthI, 0);
		if(renderParam->material.transparent)
		{
			//The dIntervalCompact is not changed.
			unfinishedN = scanDepth;			
		}
		else
		{
			//scanDepth save the number of the pixels could get the hit in further procedure.
			CUDA_SAFE_CALL( cudaMemset(dScanDepth, 0, sizeof(uint)) );
			CUDA_SAFE_CALL( cudaMemset(dIntervalCompact, 0, pixelN*sizeof(uint)) );
			CompactUnfinished<<<grid, PIXELTHREADS, 1024>>>(dIntervalCompact, dScanDepth, dUnfinished, dIntervalPad);
			CUDA_SAFE_CALL(cudaMemcpy(&unfinishedN, dScanDepth, sizeof(uint), cudaMemcpyDeviceToHost) );		
#if INFO_OUTPUT
			printf("unfinishedN is: %u\n", unfinishedN);
#endif
			if(unfinishedN==0)	break;
		}
	}
	//printf("____________The shadow test________________\n");



	//after Shadow 3DDDA, ray/surface intersection
	Db = dim3(16,16);
	Dg = dim3( width/Db.x+1, height/Db.y+1 );
	if(width %Db.x==0)	Dg.x-=1;
	if(height%Db.y==0)  Dg.y-=1;	
	Interval3DDDA2<<<Dg, Db>>>(asi, dRayDir, dDepth3DDDA, dInterval, maxInterval, dPatchIndex, nPatch, pixelBits, patchBits, (unsigned char*)surface, pitch, (int)onlyClimb, groundPlane, 1, surfNorm);

	depthI = 0;
	unfinishedN = pixelN;
	sortPadNum=0;
	grid = dim3(1,1,1);

	CUDA_SAFE_CALL( cudaMemset(dIntervalCompact, 0, pixelN*sizeof(uint)) );
	while(true&&depthI<renderParam->nIter)
	{
		grid.x = unfinishedN/SCAN_THREADN+1;
		if(unfinishedN%SCAN_THREADN==0)	grid.x-=1;
		CUDA_SAFE_CALL( cudaMemset(dScanDepth, 0, sizeof(uint)) );
		PeelDepth<<<grid, SCAN_THREADN, 1024>>>(unfinishedN, dDepth3DDDA, dInterval, maxInterval, dIntervalCompact, dScanDepth, depthI, pixelBits, patchBits);
		depthI++;

		CUDA_SAFE_CALL(cudaMemcpy(&scanDepth, dScanDepth, sizeof(uint), cudaMemcpyDeviceToHost) );
#if INFO_OUTPUT
		printf("%u step peel %u pixels.\n", depthI, scanDepth);
#endif	

		if(scanDepth==0)	break;	

		//Pad the interval to PIXELTHREADS-aligned
		dIntervalPad = dIntervalPeel;
		sortPadAll(scanDepth, &sortPadNum,  pixelN*3/2*sizeof(uint), dIntervalCompact, dIntervalPad, dIntervalPad, hIntervalSort, hIntervalPad, PIXELTHREADS);
#if INFO_OUTPUT
		printf("actual pad: %u, pad allocate: %u\n", sortPadNum, pixelN*3/2);
#endif
		//dumpBuffer(dIntervalPad, sortPadNum);

		if(sortPadNum==0)	break;
		if((int)(sortPadNum)>pixelN*maxAvgInterval)
		{
			printf("sortPad number > interval compact!\n");
			return;
		}

		//Hit test
		grid.x = sortPadNum / PIXELTHREADS + 1;
		if(sortPadNum%PIXELTHREADS == 0)		grid.x -= 1;
		CUDA_SAFE_CALL( cudaMemset(dUnfinished, 0, pixelN*3/2*sizeof(uint)) );
		PeelSolver<<<grid, PIXELTHREADS, 1024>>>((unsigned char*)surfNorm, (unsigned char*)surface, pitch, degSum, dIntervalPad, dRayDir, sortPadNum, asi, pixelBits, patchBits, dPatchIndexInverse, dFuncData, dUnfinished, depthI);
		if(renderParam->material.transparent)
		{
			//The dIntervalCompact is not changed.
			unfinishedN = scanDepth;			
		}
		else
		{
			//scanDepth save the number of the pixels could get the hit in further procedure.
			CUDA_SAFE_CALL( cudaMemset(dScanDepth, 0, sizeof(uint)) );
			CUDA_SAFE_CALL( cudaMemset(dIntervalCompact, 0, pixelN*sizeof(uint)) );
			CompactUnfinished<<<grid, PIXELTHREADS, 1024>>>(dIntervalCompact, dScanDepth, dUnfinished, dIntervalPad);
			CUDA_SAFE_CALL(cudaMemcpy(&unfinishedN, dScanDepth, sizeof(uint), cudaMemcpyDeviceToHost) );		
#if INFO_OUTPUT
			printf("unfinishedN is: %u\n", unfinishedN);
#endif
			if(unfinishedN==0)	break;
		}
	}


	//Edge detect of the render result for anti-aliasing:
	if(renderParam->renderState==1&&renderParam->material.transparent==0)
	{
		//printf("----------------------------\n");
		uint* dAAPointsN = dScanDepth;
		uint AAPointsN = 0;
		uint* dAAFlag = dDepth3DDDAPeel;
		uint* dAAPoints = dIntervalCompact;
		uint* dRSAAPoints = dIntervalPeel;

		CUDA_SAFE_CALL(cudaMemset(dAAPointsN,  0, sizeof(uint)));
		CUDA_SAFE_CALL(cudaMemset(dAAFlag,     0, pixelN*sizeof(uint)));
		CUDA_SAFE_CALL(cudaMemset(dAAPoints,   0, pixelN*sizeof(uint)));
		CUDA_SAFE_CALL(cudaMemset(dRSAAPoints, 0, pixelN*sizeof(uint)));

		SilhouetteEdgeDetect<<<Dg, Db, 1024>>>((unsigned char*)surfNorm, (unsigned char*)surface, pitch, dAAFlag, dAAPointsN, asi, dPatchIndex, nPatch, pixelBits, patchBits);

		Shading<<<Dg, Db>>>((unsigned char*)surfNorm, (unsigned char*)surface, pitch, asi);

		CUDA_SAFE_CALL(cudaMemcpy(&AAPointsN, dAAPointsN, sizeof(uint), cudaMemcpyDeviceToHost) );
		//dumpBuffer(dAAFlag, AAPointsN);

		sortPadAll(AAPointsN, &sortPadNum, pixelN*sizeof(uint),      dAAFlag, dRSAAPoints, dRSAAPoints, hIntervalSort, hIntervalPad, PIXELTHREADS/RSAA_SAMPLE);
		//dumpBuffer(dRSAAPoints, sortPadNum);
		if(sortPadNum==0)	return;

		sortPadNum *= RSAA_SAMPLE;
		grid.x = sortPadNum / PIXELTHREADS + 1;
		if(sortPadNum%PIXELTHREADS == 0)		grid.x -= 1;
		RSAARayCasting<<<grid, PIXELTHREADS, 1024>>>((unsigned char*)surface, (unsigned char*)surfNorm, pitch, degSum, dRSAAPoints, dRayDir, sortPadNum, asi, pixelBits, patchBits, dPatchIndexInverse, dFuncData);
	}

#if SELF_SHADOW
	//Self Shadow
	Db = dim3(16,16);
	Dg = dim3( width/Db.x+1, height/Db.y+1 );
	if(width %Db.x==0)	Dg.x-=1;
	if(height%Db.y==0)  Dg.y-=1;	
	Interval3DDDA2<<<Dg, Db>>>(asi, dRayDir, dDepth3DDDA, dInterval, maxInterval, dPatchIndex, nPatch, pixelBits, patchBits, (unsigned char*)surface, pitch, (int)onlyClimb, groundPlane, 2, surfNorm);

	depthI = 0;
	unfinishedN = pixelN;
	sortPadNum=0;
	grid = dim3(1,1,1);

	CUDA_SAFE_CALL( cudaMemset(dIntervalCompact, 0, pixelN*sizeof(uint)) );
#if INFO_OUTPUT
	printf("Self shadow begin.............................\n");
#endif
	while(true&&depthI<renderParam->nIter)
	{
		grid.x = unfinishedN/SCAN_THREADN+1;
		if(unfinishedN%SCAN_THREADN==0)	grid.x-=1;
		CUDA_SAFE_CALL( cudaMemset(dScanDepth, 0, sizeof(uint)) );
		PeelDepth<<<grid, SCAN_THREADN, 1024>>>(unfinishedN, dDepth3DDDA, dInterval, maxInterval, dIntervalCompact, dScanDepth, depthI, pixelBits, patchBits);
		depthI++;

		CUDA_SAFE_CALL(cudaMemcpy(&scanDepth, dScanDepth, sizeof(uint), cudaMemcpyDeviceToHost) );
#if INFO_OUTPUT
		printf("%u step peel %u pixels.\n", depthI, scanDepth);
#endif	

		if(scanDepth==0)	break;	

		//Pad the interval to PIXELTHREADS-aligned
		dIntervalPad = dIntervalPeel;
		sortPadAll(scanDepth, &sortPadNum,  pixelN*3/2*sizeof(uint), dIntervalCompact, dIntervalPad, dIntervalPad, hIntervalSort, hIntervalPad, PIXELTHREADS);
#if INFO_OUTPUT
		printf("actual pad: %u, pad allocate: %u\n", sortPadNum, pixelN*3/2);
#endif
		//dumpBuffer(dIntervalPad, sortPadNum);

		if(sortPadNum==0)	break;
		if((int)(sortPadNum)>pixelN*maxAvgInterval)
		{
			printf("sortPad number > interval compact!\n");
			return;
		}

		//Hit test
		grid.x = sortPadNum / PIXELTHREADS + 1;
		if(sortPadNum%PIXELTHREADS == 0)		grid.x -= 1;
		CUDA_SAFE_CALL( cudaMemset(dUnfinished, 0, pixelN*3/2*sizeof(uint)) );
		PeelSolver2<<<grid, PIXELTHREADS, 1024>>>((unsigned char*)surfNorm, (unsigned char*)surface, pitch, degSum, dIntervalPad, dRayDir, sortPadNum, asi, pixelBits, patchBits, dPatchIndexInverse, dFuncData, dUnfinished, depthI, 2);
		if(renderParam->material.transparent)
		{
			//The dIntervalCompact is not changed.
			unfinishedN = scanDepth;			
		}
		else
		{
			//scanDepth save the number of the pixels could get the hit in further procedure.
			CUDA_SAFE_CALL( cudaMemset(dScanDepth, 0, sizeof(uint)) );
			CUDA_SAFE_CALL( cudaMemset(dIntervalCompact, 0, pixelN*sizeof(uint)) );
			CompactUnfinished<<<grid, PIXELTHREADS, 1024>>>(dIntervalCompact, dScanDepth, dUnfinished, dIntervalPad);
			CUDA_SAFE_CALL(cudaMemcpy(&unfinishedN, dScanDepth, sizeof(uint), cudaMemcpyDeviceToHost) );		
#if INFO_OUTPUT
			printf("unfinishedN is: %u\n", unfinishedN);
#endif
			if(unfinishedN==0)	break;
		}
	}
#if INFO_OUTPUT
	printf("Self shadow end.............................\n");
#endif
	printf("The anti-aliasing of self-shadow is not handeled!\n");
	printf("We need receiver-Shader + ray casting + self-Shader, and anti-aliasing in the last step!\n");
#endif
}